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(March 25, 1964) 

An iterative procedure for unfolding the effects of ttie finite resolution of a detector 
from an observed pulse height distribution is discussed. The process is demonstrated for a 
particular detection system. Convergence and uniqueness properties of the method are 
discussed empirically. 

A general expression for the propagated error resulting from errors in the detected 
pulse height distribution is derived. Approximations are made in order to evaluate the 
propagated error for a particular detector. These approximations become better as the 
resolution of the detector Improves. The results indicate that the error rapidly approaches 
a limit of from 1.5 to 3 times the error in the observed distribution. This limit is reached 
in approximately three iterations. 



1. Introduction 

In measiiriiio' y-i'iiy spectra it is frequently neces- 
sary to remove tlie effects of the resolution of tlie 
detector from an observed pidse height distribution. 
This is known as ^' unfolding'^, '^mscrambling'', or 
^'unsmearing^'. To do this a nuitrix representing 
the response of the detector must be found. Let 
the incident spectrum be denoted by an rn dimen- 
sional vector A^: 






N= 



The response may be represented by an mXw 
matrix R. The detected pulse height distribution, 
P, is then given by 



Pi=i: iiuNi. 



(1) 



Unfolding is the name given to the process of finding 
A^ such that 

j 

where Rij~^ is the inverse to the matrix Rtj. 

It is frequently undesirable to obtain a solution 
Nj by inverting the response function matrix. 
Usually the response function matrix is a very large 
square matrix. In this experiment one form of the 
matrix was 700X700. The inversion of such a 
matrix would be a formidable task, even when 
utilizing computer techniques. 

For this reason, iterative approximations to solu- 
tions have been developed by Scofield |2J and hy 
Skarsgard, Johns, and Green [:^J. An iterative 
technique similar to that described by tlie latter 
has been developed inde[)eiulently in this laboratory. 
Conver2:ence criteria lor (his technique have been 
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discussed by Geiringer [1]. In applying the tech- 
nique, empirical evidence has been obtained for the 
vahdity of solutions obtained by this method. This 
evidence is discussed below. In addition, the 
propagation of error in the unfolding process is 
investigated in detail. 

2. Iterative Solution 

2.1. Procedure 

Equation (1) may be written in matrix form as 

P=RN. 

Assume N= T/y. Then this initial estimate will give 

A measure of the closeness with which U^ represents 
the true A^ is given by the difference 

Ao-P-Po. 
Uq may be corrected to form 

and the new correction 

Ai-P-PC/i 
is found. For the in}^ iteration 

Un^Un-^ + ^n-^ 1 

A.> = F-/nJ, . (2) 

It has been found for the present work that it is 
satisfactory to use P itself as the initial estimate Uq. 
The technique has been used primarily in unfolding 
pulse height distributions obtained witli the NBS 
Two Crystal Pair Spectrometer. Details of the 
detector and its response are described by Ziegler, 
Wyckoff, and Koch [4]. 

Various methods for arresting the iterative pro- 
cedure may be used. In this work the data were 
unfolded using a predetermined number of iterations. 

2.2. The Response Function Matrix 

This section will (hscuss the problem of finding 
a matrix representation for the assumed analytic 
form of the response function. The response at 
pulse height e due to one incident photon of energy 
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k^ may be written [4] 



(3) 



J E{e,h)dko=] 



where Ci, C2, C^, d, and a are constants character- 
istic of the detector. The pulse height distribution 
becomes 



^«=r 



dhR(e,ko)N(h) 



(4) 



where N(ko) is the continuous incident photon num- 
ber spectrum. Equation (4) is the continuous form 
of (1). Experimentally the vector 

Pi= P(e)de 

is the quantity measured as counts per channel in 
a multichannel pulse height analyzer. 

The integral equation (4) does not possess an 
exact solution. Integrating over k in (3), (4) 
becomes 

PW =jdko I aC, exp [-^^5^'] 

+a^.("Vl)[-p(-^+zi 



+ 






where ^(x) is the error integral [5]. Thus (5) is 
seen to be an integral equation with a Gaussian 
kernel; such an equation does not possess a general 
unique solution [6]. This is a manifestation of the 
inability to experimentally differentiate between a 
smooth spectrum and a spectrum containing a series 
of sharp spikes. The Gaussian broadening is respon- 
sible for this. 

In order to obtain a matrix representation for 
R(€, ko) a particular form must be assumed for 
Niko). Two forms have been investigated. One 
may assume the incident spectrum to consist of a 
series of discrete steps so that over a fixed small 
energy width the spectrum is constant [7]. Alter- 
natively one may assume the spectrum to be com- 
posed of a sum of Dirac delta functions so that 
when an integration is performed over a small energy 
width the area of the delta function gives the number 
of photons in that width. 



Both cases lead to essentially the same form for 
the matrix. The latter case will be carried through 
to obtain the matrix explicitly. 



Let N(ko)=^aj8{kj—ko). Then (4) becomes 
P(e)=j:Kjm,e)aj 



(6) 



where Kj=CiCA+kjC2C3Ci(l — e ^^^3) is a number, and 

J{kj,e)=e-<—)' 

From the above the number of counts in channel 

e^is 



Pi= r^^' J:Kjf(kj, e)aj 

Jei-Ai j 



which, after interchanging integration and sum- 
mation becomes 



Pi=j:aj[Kjb(kj,e,,A,)] 



(7) 



where 



b{kj,e„Ad=\ e ^- - Ue = bj,(A^). 

J ei-Ai 

Identifying aj with Nj and Kjbu with Rfj (7) becomes 
identical with (1). 

3. Empirical Justification 

3.1. Convergence 

In setting out on this course there was no reason 
to believe the technique to be convergent. It has 
been shown [3] that convergence is assured for a 
smooth function if the eigenvalues A^ of the response 
function matrix satisfy the requirement 

0<A,<2. 

This was not a useful test because the size of the 
matrices used made calculation of the eigenvalues 
impractical. Therefore the primary justification 
is empirical. 

In analysis utilizing a 200 X 200 form of Rij, eleven 
iterations were ordinarily performed. However, 
as a check on convergence, as many as twenty-one 
iterations have been performed, during which A^ of 
equation (2) continues to converge. 

In figure 1 a typical set of points to be unfolded is 
plotted. Let A denote this set. On the same figure 
is plotted B, the result of unfolding A. The set A 
contained points only up to 40 MeV. In order to 
avoid introducing a large discontinuity in the first 
derivative at 40 MeV, a straight line tail has been 
added. The work was done with an energy grid 
width of 0.5 MeV. Typical standard deviation is 
shown for a point of ^ at 16 MeV. 

In order to compare B with A, the difference 
A^^=A-RB (see (2)) is plotted in figure 2. If B is 
the ^'correct'' unfolded set of points then An must 
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Figure 1. Original curve^and the^corres'ponding unfolded 
curve after eleven iterations, for ji typical spectrum. 
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Figure 2. Difference, A, between the original curve of figure 1 
and the folding of the unfolded curve of figure 1, 

X% indicates percentage dilTcrencc A/.l. Ordinate scale is in same units as 
figure 1. 



vanisli. ronvergenco roqiiii-es that z^,^ vanishes for 
increasing n. 

Some values of (lie (HO'erence in percent are in- 
dicated on the i^lot. The very small (0.7%) differ- 
ence at 19.5 MeV is at the peak of A. 

In order to furtlier check the convergence proper- 
ties of the scheme a set of points with large 
uncertainties was unfolded, using twenty-one itera- 
tions. The set C and its 'Hmfokr^ D are shown in 
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Rt 3. Original curve with poor statistics and unfolded 
curve after twenty-one iterations. 
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Figure 4. The difference An between the original curve of 
figure 3 and the folding of the unfolded curve after n iterations, 
for n=0, ly and 4- 

Ordinate scale in figures 4 and 5 is in the same units as figure 3. 

figure 3. Again a straight line tail has been added 
to C at 40 MeV. 

Because of tlie poorer statistics on C there are 
more fluctuations in D. The slope of C appears to 
have a large discontinuity at 32.5 MeV. A spike in 
D is observed to grow at this energy with successive 
iterations. This demonstrates that fluctuations are 
magnified as one approaches an exact solution. 

The question of convergence is best illustrated by 
examining A^^, for various values of n. Figures 4 
and 5 show A^ for 71=0,1,4,11, and 21. It is ob- 
sei'ved that A„ con\erges rapidly for small ii. The 
maximum of the ratio AoJAq is approximately 10"^. 
Tlie maxima of A21 in percent of C are +0.4 and —0.5 
[)ercent. 

A numerical criterion for testing convergence in 
this sense is suggested by Skarsgard, Johns, and 
Green [3] for a pulse height distribution containing 
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Figure 5. The difference Aq for n=4, H, and 21. 
Note that the ordinate scale is expanded from figure 4. 



pure Poisson counting errors, i.e. the standard 
deviation on P^ is V^r 

If p <^<^1 then the deviation for the point 

* . . 

Ui is well within the limits of random measm^ement 

errors. Therefore if 



M a2 
i = l i^i 



(8) 



the unfold is regarded as satisfactory. This test 
was used in unfolding a pulse height distribution 
for which the errors on each point were purely count- 
ing errors. The results were similar to those found 
by Skarsgard, Johns, and Green, [3] namely, con- 
vergence is rapid until (8) is satisfied (^3 itera- 
tions). After this convergence proceeds slowly. 

One might hope to be able to prove convergence 
from the classical theorems [1]. It is easily shown 
that if one denotes {I-R) by A, then: 

Therefore [1] ?7^"^ converges to the true solution 
Ut if and only if the eigenvalues of A are less than 
one in modulus. From the rapid convergence which 
is observed empirically, one is led to believe that 
the eigenvalues of A are indeed less than one in 
modulus. 

A sufficient condition for convergence is that 
the maximum of the absolute row sums ai, satisfy 



M.=(|]|A,|) 



<1. 



However, this is not the case for the matrix A=I-R 
on which the present work is based. 



4. Error Propagation 
4.1. Empirical 

In order to demonstrate the effect of statistical 
fluctuations, two different experimental determina- 
tions of the same pulse height distribution have 
been unfolded. A portion of the unfolded spectra 
for both sets of data are presented in figure 6. 
The spectra are designated Ua and Ut,. The two 
pulse height distributions are not presented because 
of the typographical difficulty in distinguishing the 
two sets of data on a meaningful scale. 

Differences between the two spectra should be 
purely statistical. Let the measured pulse height 
distributions from which Ua and U^ were obtained 
be Pa and P,. The ratio p={Pa/Ua)KP,IU,) has 
been plotted in figure 7 for the region from 15 MeV 
to 25 MeV. One would expect this ratio to be ran- 
domly distributed about imity due to the statistical 
fluctuations in Pa and P^. This is observed. 

In addition, if the unfolding procedure does not 
introduce false structure, then, for Pa^Ph one ex- 
pects the relations Ua^Uj^ and therefore p>l, to 
hold approximately. Examination of figures 6 and 
7 will show that for Ua^Uj,, p>l, and for Ua<^Uf), 
p<^l, except when Ua^^^Uf, where fluctuations in 
adjacent points become important. This indicates 
that the iterative procedure does not introduce false 
structure. 

Some qualitative eflects of error propagation are 
illustrated quite well by figure 6. The most pro- 
nounced effect is the increase of fluctuations in the 
unfolded curves with increasing energies. The rea- 
son for this will emeroe from the discussion followino". 
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Figure 6. A portion of the unfolded spectra for two different 
experimental determinations of the same pulse height distri- 
bution after eleven iterations. 
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Figure 7. The ratio p= (PJV,)/(FJ\J^). 
See text for definitions. 

4.2. Calculation of Error 

The calculation of error for an individual point in 
an unfolded spectrum is made difficult, because, in 
folding, correlations arise between errors in adjacent 
points. 

Assume that error in the detected pulse height 
distribution is known as a function of energy and 
denote it by 0-^- The folded set of points, which are 
obtained in the first step, may be written 



The standard deviation of Sf becomes 



(9) 



(10) 



From (2) the solution after one iteration may be 
written 



Ui' 



--2Pi-'EIiijPj 



with corresponding error 
The second iteration gives: 



(11) 



(12) 



where 



=SP-3EP+R'P 
R^P=^RijRjkPk' 

jk 



Expanding U^^^ in the same way as U^^\ the error on 
U^^^ becomes 

8Ui^'^ = [(S-3Ru+Ru)<Ti'+ other terms]^. 

The ^'other terms^' have numerous cross products. 



For example, the contribution to [5?7P]^ from 0-^4-1 is 

[ ^Ri,i+l~{~RiiRi,i+l~{'Ri,i + lRi+Ui+l] O'i-f-l • 

The general form for the solution after n iterations 
may be written symbolically as 



C/^-^=ir/_(/_i?)n+i1p 



(13) 



where is the identity matrix and 1/R=R~\ If 
the variance Var (P) = (t^I, then the variance of 
?7(^) may be written formally as: 

Var (f/^^O 

= R-'[I-(I-Rr+'][I-(I-Ry-^T(R~Yc^\ 

where T denotes the transpose. Here we have used 

[8] 

Var(f7P) = (7[var (P)]C'^. 

Since all elements of 72 are less than unity, it is evi- 
dent that in the limit as n approaches infinity 



Um 



Var(f7^-0=^"'(^"0V. 



In order to simplify the error calculation the 
response function and the error will be assmned to 
satisfy the following conditions: 

(1) The half -width of tlie response function is 
narrow. This corresponds to good resolution in the 
detector. 

(2) The shape of the response function does not 
change rapidly with incident photon energy. This 
is equivalent to assuming tliat 

R(e,h)^Rie+Akuki + Ah) 



or Rij^Ri^ 



and Ri_.rn,j^Ri.3 



where m is an integer. See figure 8. 

(3) The error, c^., is a constant, a, over the half- 
width of the response function. Note that the first 
condition makes this more likely. 

het Rii=Wo, Ri-i,i=Wij Ri-2,i='W2, • • • But, from 
condition (2) Ri,i^j,o^Wic. Using this (9) may now be 
written 



[RP]i=Jlw,P,+i, 



(9a) 



where the w^, may be obtained from the response at 
one incident energy. 

In the appendix it is shown that the solution has 
the o:eneral form: 



Ul-'=a,P,+a,P,^,+a2Pi+2+ 



(14) 



If the error is assumed to be a constant, a, this 
arives 






(15) 
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Figure 8. Illustration of the relation between the response 
functions at two different energies. 



The folded spectrum at i is RuPi-^JRi ,+iP,+i-|- . . 
in the text. Bi i+i^Ri- 



. where R and P are defined 
■i.i. 



for the error at point i in the nth iterated solution. 

For the response described in [4] the error after 
three iterations has been calcuhited at 18.5 MeV and 
48.5 MeV. (See appendix.) The results are ±2.36 
a at 18.5 MeV and ±3.23 a at 48.5 MeV. 

In both cases if all terms in (A3) which are cubic 
in w were omitted, the difference in C7/^^ would be 
small, and the difference in 517/^^ would be negligible. 
The terms which are cubic in w are approximately 
an order of magnitude smaller than the quadratic 
terms. The conclusion is that §?7^"^ has converged 
forn>3. 

From (^1), (A2), and {A3) it may be observed 
that after a large number of iterations tlie coefficient 
of Po in {A3) will converge to: 



tto 



4l-(l-^o)^+^]M. 



For the cases at 18.5 and 48.5 MeV this gives 
ao=2.21 and 3.104 respectively, for n=3. Compari- 
son with the results for 5f7/^^ above shows very close 
agreement. Thus one concludes that three iterations 
satisfy the large number criterion. 

This coefficient then places a lower limit on the 
propagated error at each point of the solution. In 
general Wo decreases with increasing incident photon 
energy, for this experiment. Therefore the error 
must increase with energy. This is independent of 
the shape of the curve to be unfolded. Figure 9 
shows Go plotted as a function of energy for n=3. 
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Figure 9. The function ao = [l — (1 — Wo) °+^]/wo as a function 
of energy for n=3. 



5. Summary 

Because of the uncertainties associated with any 
point on a measured pulse height distribution, any 
' 'solution'' for the unfolded spectrum is acceptable, if 
the difference between the measured distribution and 
the fold of the solution lies within the uncertainty 
associated with the measured distribution. The 
additional requirement of smoothness is sufficient to 
ensure that the iterative process converges to a 
useful solution. 

General error analysis is difficult. However, 
approximations may be made which become better 
as the resolution of the detector improves; these 
approximations make an error estimate possible. 

Interesting results from the use of this technique 
may be seen in tlie work of Ziegler, Koch, Wyckoff, 
and Uhlig [9]. 
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encouragement of H. W. Koch, R. A. Schrack, and 
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and their suggestions for improvement of this paper. 
The original work in finding a suitable form for the 
response function and the fitting of that form to 
experimental data was done by R. A. Schrack. 

Special thanks go to Joseph Cameron and Brian 
Joiner for their help in studying the propagation of 
error. 



6. Appendix 

In this appendix it will be shown how the general 
form of the solution, equation (14), may be found. 
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Using (9a) and expanding (13) the general solution 
becomes: 



2! 

nin—l)(:n-2) 



3! 



k j 



7i{7h—l){n-2)('n-3) 
4! 

J2^^^i^\'^^jPi+fc+j+i+ . . . (Al) 



Choosing an arbitral'}'' zero point for index i (Al) 
may be written: 

- J^w^P.j - 3 2 Z;^I^,-P2.+ - 3 S ^wlwjP2,-,j 

j k<j j<k 



~~G S XI 2^'^^+^•+.^ 
I <k<j 



(A2) 



Expanding and cc^llecting terms: 

U^^^ = {4—6iOo+Awl—wl)Pi^-{- (—6wi+HwoWi—'3woWi)Pi 
+ ( — 6w2 + 4wl + HW0W2 — 3wlw2 — 3WiiWl) P2 
-{- (—(jw-s+SwoW3+SwiW2—wl—3wlw3 
— 6WQW1W2) Pi + (—(JW4 + 4:W2 + 8i(;o'W^4 + 8'w;it(;3 
— dwiwi — l^wlw2 — SwIwq — QW0W1W3) P4 + ( — 6ic'5 
+ SwqWs + 8t(;ii(;4 + Sw2Ws — ^wlws — ^ivlw^ — Swlwi 
—QwoWiWi—QwoW2Ws)P5-{- . . , (A3) 

At 18.5 MeV the response function for the detec- 
tion system described by Ziegler, Wyckoff, and 
Kocli [4] was given by 



(Wo,Wi,W2,Ws,W4,Ws, . . .) 

= (0.39,0.225,0.138,0.088, 

0.056,0.035, . . 

Using these values in (A3) one finds: 

C7o ^'^=2.209 7^0-0.751 Pi-0.317 P2 

-0.129 P3-O.O84 P4-O.OO3 P5- . 



and from (15): 



dW^ =2.^6(7. 



Similarly the response at 48.5 MeV gives 

Z7o^^) = 3.103 Po-0.642 Pi-0.454 P2 

-0.333 P3-O.242 P4-O.174 7^5- 



which leads to 



6[7o^'^=3.23(r. 
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